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Single-particle resonant states, also called Gamow states, as well as bound and scattering states 
of complex energy form a complete set, the Berggren completeness relation. It is the building block 
of the recently introduced Gamow Shell Model, where weakly bound and resonant nuclear wave 
functions are expanded with a many-body basis of Slater Determinants generated by this set of 
single-particle states. However, Gamow states have never been studied in the context of Hartree- 
Fock-Bogoliubov theory, except in the Bardeen-Cooper-Schriefer (BCS) approximation, where both 
the upper and lower components of a quasiparticle wave function are assumed to possess the same 
radial dependence with that of a Gamow state associated with the Hartree-Fock potential. Hence, 
an extension of the notion of Gamow state has to be effected in the domain of quasiparticles. It 
is shown theoretically and numerically that bound, resonant and scattering quasiparticles are well 
denned and form a complete set, by which bound Hartree-Fock-Bogoliubov ground states can be 
constructed. It is also shown that the Gamow-Hartree-Fock single-particle basis can be used to solve 
the Gamow-Hartree-Fock-Bogoliubov problem. As an illustration, the proposed method is applied 
to neutron-rich Nickel isotopes close to the neutron drip-line. 
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I. INTRODUCTION 



One of current challenges of nuclear theory is the quantitative description of nuclei situated near and beyond drip- 
lines. Powerful facilities are being built in several countries in order to create these very short-lived states. For a long 
time, microscopic theories of nuclear structure have been developed for describing ground states of nuclei close to the 
valley of stability. For describing stable nuclei which are well localized, the harmonic oscillator (HO) bases are useful 
for both shell model [l[ and Hartree-Fock Bogoliubov (HFB) calculations 0,0, 0, [H, [|| ; the HO bases converge quickly 
therein. However, it possesses poor convergence properties for weakly bound nuclei bearing large spatial extensions, 
which lie very close to neutron drip lines. 

A promising approach to this problem has been proposed in Refs. @, HHl within a shell model framework; namely, 
the Gamow Shell Model (GSM). The fundamental idea is to replace the HO basis by the Berggren basis consisting 
of bound states, resonance states and continuum scattering states of complex energy, generated by a single-particle 
potential. It has been shown numerically that this basis has the ability to expand both halo nuclei and many-body 
resonant states precisely. The latter indeed belongs to a rigged Hilbert space [ToL [Tlj , which is an extension of the 
notion of Hilbert space to non-square integrable wave functions. However, the dimension of the Berggren Slater 
Determinants represented by the GSM basis increase very quickly with increasing number of valence particles; it 
increases much faster than in standard shell model due to the presence of occupied scattering states. Hence, the GSM 
is a tool mainly dedicated to the study of light nuclei. For medium and heavy nuclei, a method of choice is the HFB, 
which can be followed by quasiparticle random phase approximation (QRPA). As pairing correlations are absorbed in 
the HFB ground state, one-body nature of the HFB framework enables fast evaluations of ground states of medium 
and heavy nuclei, and it is in fact the only method suitable for systematic calculations; see Ref. [l2j for an evaluation 
of even-even nuclei in the whole nuclear chart with the HFB formalism. In order to properly treat drip-linc nuclei 
within the HFB framework, the real-space coordinate-mesh method has been applied using box boundary conditions 
[HI, [HI ■ Extension of this approach to deformed nuclei is difficult and has been carried out only recently [HI, [n| • As 
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an alternative more convenient approach, one can adopt basis expansion methods, where direct integration procedure 
is replaced by matrix diagonalization. A first amelioration of the HO basis had been proposed with the transformed 
harmonic oscillator (THO) basis [H, E3|- Applying unitary transformations to the HO basis, one obtains the THO 
basis, in which Gaussian fall-off of the HO wave functions is replaced by physical exponential decrease of the THO 
basis wave functions. However, the THO basis always dictates exponential decrease in expanding quasiparticlc wave 
functions, for both upper and lower components, even when they are part of scattering states, so that unsatisfactory 
basis dependence remains. In order to solve this problem, a new basis has been introduced very recently, which consists 
in using bound and continuum basis statesgenerated by the analytic Poschl-Teller-Ginocchio (PTG) potential [l8l |. 
The PTG basis introduced in this paper [19[ possesses a peculiarity to bear no narrow resonance states; those are 
replaced by bound PTG states. Thus, PTG continuum set of basis states can be discretized very effectively with 
Gauss-Legendre quadrature, as they contain no resonant structure. It has been shown that they can provide a good 
description of spatially extended nuclear ground states of both spherical and axially deformed nuclei On the 

other hand, the PTG basis formed by bound and real scattering states is not a Berggren complete set of states, so 
that it would be more convenient to use a Berggren quasiparticle basis set, when we are interested in describing 
particle-decaying excited states. Up to now, however, resonant quasiparticle states have been studied in the context 
of Berggren completeness relation only within the BCS approximation [2^, l2lj . The last approach is indeed not 
satisfactory due to the well-known gas problem arising from the occupation of the continuum: In fact, densities are 
not localized in the BCS approach, because the lower components of scattering quasiparticle states are of scattering 
type as well. Contrary to what is stated in Ref. [201 ] . it cannot be regularized using complex scaling because it does 
not have pure outgoing asymptotic. Use of continuum level density in Ref. [2 if is also problematic, even though it 
suppresses the gas problem. Indeed, it is not part of continuum HFB theory [lj|, so that its introduction in HFB 
equations strongly modifies quasiparticle coupling to the continuum. In particular, it supresses a large part of non- 
resonant continuum, and thus important physical properties of drip-line nuclei as well. Hence, with this approach, 
weakly bound systems cannot be studied properly. Only a full application of the HFB framework can unambiguously 
solve this problem, where densities are localized by construction for bound HFB ground states. 

The major purpose of this paper is to develop a new method of solving the continuum HFB equations utilizing 
the Berggren basis, called Gamow-HFB method, by which bound, resonant and continuum quasiparticle states are 
provided. It allows expansion of QRPA excited states having escaping widths in terms of the Berggren quasiparticle 
basis associated with the bound HFB ground state. This is very important because, in weakly bound unstable nuclei, 
low-lying collective excited states may acquire particle-decay widths. 

This paper is organized as follows. Firstly, the standard HFB formalism is briefly summarized. As we use the Skyrme 
interactions [13], it is effected in the context of density functional theory (DFT). Secondly, we define quasiparticle 
S-matrix poles and scattering states of complex energy; these arc direct extensions of their single-particle counterparts. 
We then present the quasiparticle Berggren completeness relation generated by those states. Numerical methods to 
calculate Gamow and complex scattering quasiparticle states are described; they differ significantly from the scattering 
quasiparticle states discretized by box boundary conditions. We also present another method of solving the continuum 
HFB equations in which the HFB quasiparticle wave functions are expanded in terms of the Gamow-Hartree-Fock 
(GHF) basis; this approach may be regarded as an extension of the standard two-basis method [23|,[24], to complex 
energy plane. Feasibility of the proposed methods is illustrated for neutron-rich Nickel isotopes close to the drip line. 
Perspectives for unbound HFB theory and QRPA calculations using the Gamow-HFB quasiparticle basis will then be 
discussed. 



II. GENERAL HFB FORMALISM WITH DFT 



The HFB equations are expressed in super-matrix form constituted by particle-hole field Hamiltonian h, particle- 
particle pairing Hamiltonian h and chemical potential A guaranteeing conservation of particle number in average: 

Y .-„)(:)=<:)• w 

Using Skyrme and density-dependent contact interactions for the particle-hole and pairing channels, respectively, h 
and h are expressed in terms of local normal density p{r) and pairing density p{r). Formulas providing p, p, h and h 
can be found in [flU^I- As h and h depend on p and p, determined from quasiparticles eigenvectors of Eq. (fTJ), the 
HFB equations must be solved in a self-consistent manner [27j . 

Let us consider the HFB equations with the Skyrme energy density functionals and density-dependent contact 
pairing interactions assuming spherical symmetry. Fixing orbital and total angular momentum £ and j, as well as 
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proton or neutron nature of the wave functions, Eq. ((T|) becomes a system of radial differential equations [lj 



dr 2m* (r) dr 
d h 2 d 
dr 2m* (r) dr 



u(k, r) = 
v(k, r) = 



h 2 £(£ + l) 
2m*(r)r 2 

n 2 £(e + i) 

2m*(r)r 2 



+ V(r) - (A + E) 
+ V(r) - (A - E) 



u(k, r) + W(r) v(k, r), 
v{k, r) — W(r) u(k, r), 



(2) 



where 



• u(k,r) and v(k,r) are respectively the upper and lower components of quasiparticle wave function with energy 
E, and k = \/2mE '/h with the nucleon mass m, 

• to* (r), V(r) and W(r) are respectively the effective mass, the particle-hole (field) and particle-particle (pairing) 
potentials of the HFB Hamiltonian. 



Because nuclear interactions are finite range, only Coulomb and centrifugal parts remain for r 
becomes asymptotically: 



-co, so that Eq. © 



d 2 u 



dr 2 
dr 



(fc,r) 



{1(1 + 1) , 2 Vu k u 



2 {k,r) 



( £(t + l) 

I r 2 



+ 



r 

2rj v k v 



u(k,r), 
v(k,r), 

where the generalized momenta k u , k v and their associated Sommerfeld parameters 77 u , rj v are defined by 

kn, — 



2m,, „. 



mZC c 



h 2 k, 



(v) 



(proton) , rj u ^ = (neutron) 



(3) 

(4) 
(5) 



with the number of protons Z and the Coulomb constant C c . Hence, u(k,r) and v(k,r) are linear combinations of 
the Hankel or Coulomb wave functions H ir) ( ) (fc u (u) r ') for r — ► +oo. Note that k v is always imaginary provided the 

HFB ground state is bound (A < 0), while k u is real (imaginary) for E > — A (E < — A). 

The chemical potentials A for neutrons and protons are determined from the requirement of conservation of their 
number in average: 



N , Ni 



+ OC 



v i (r) dr, 



(6) 



(and similar equations for protons). Here the sum runs over all single-particle states, N is the number of neutrons 
and (N) is the expectation value in the HFB ground state. For a given particle-hole field Hamiltonian h, the chemical 
potential A could be calculated in principle exactly at each iteration, recalculating all quasiparticle wave functions from 
Eq. {T]) and updating A until Eq. (J5]) is verified. However, in practice, it is much faster to use instead an approximate 
chemical potential issued from the BCS formulas, which will converge self-consistently to the exact chemical potential 
along with the HFB Hamiltonian [l4l |. For that, one defines auxiliary single-particle energies e-i and auxiliary pairing 
gaps A by 



= A + Ei(l — 2Ni) , Aj = 2E iy /N i (l-N i ), 



(7) 



which are defined by applying the BCS type formula to the HFB quasiparticle energies the average particle number 
Ni defined in Eq. ([6]) and the chemical potential A issued from the previous iteration. While e\ and Aj correspond to 
the single-particle energy and the pairing gap in the BCS approximation, they are used here as auxiliary variables to 
solve the HFB equations. The approximate chemical potential A is obtained by solving its associated BCS equation: 



E 



t - 



A 



a/ (ei - A) 2 + A? 



2N. 



(8) 
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III. S-MATRIX POLES AND SCATTERING QUASIPARTICLE STATES 



A. Boundary conditions 



The upper and lower components, u(k,r) and v(k, r), of the quasiparticle wave function satisfy the following 
boundary conditions: 

u(k, r) ~ Cy +1 , v(k, r) ~ Cy +1 , r -f (9) 

u(k, r) ~ (M + C"^ (fc u r) , r -> +^ (10) 

«(fc,r) - (M , r -> +oo. (11) 

Eq. ([9]) is required by regularity of wave functions at r = 0. Eqs. (|10p and (jTTJ) determine the nature of quasiparticle 
state, which can be a bound, resonant (C~ — 0) or scattering (C~ ^ 0) state, and are generalizations of the boundary 
conditions defining single-particle states using the Berggren completeness relation. Eq. (jTTJ) demands outgoing wave 
function behavior of v(r) for all quasiparticle states. If its energy E is real and positive, as in the standard HFB 
approach, Eq. (|11|) is equivalent to the asymptotic condition v(k,r) — > for r — > +oo; the condition arising from 
integrability of nuclear density over all space [14j. Extension to complex energies follows from analyticity of the v(k,r) 
function in the complex fc-plane. Eq. (fTOf with C~ = then defines quasiparticle S-matrix poles, as it is equivalent 
to u(k, r) — ► for r — ► +oo for bound quasiparticle states with E < |A|, and provides resonant quasiparticle states if 
E is complex. Eq. PH|) with C~ ^ represents standard scattering quasiparticle states for real and positive E, but 
they are extended to complex energies by analyticity arguments. 



B. Normalization of quasiparticle states 



Bound HFB quasiparticle states with energy E n are normalized by: 
[u(k n ,r) 2 +v(k n ,r) 2 ] dr = 1. 

o 



(12) 



where k n — y/2mE n /h. For resonant quasiparticle states, the integral in the above equation diverges, so that this 
normalization condition cannot be used. The complex scaling method has been known as a practical means to 
normalize single-particle resonance states 28] . Convergence of integrals is obtained therein integrating up to a finite 
radius R situated in the asymptotic region, after which the interval [R : +oo[ is replaced by a complex contour defined 
by a rotation angle 8 > 0, allowing exponential decrease of the integrand. Owing to Eqs. (flT))) and (TTTj) , the same 
method can be used to normalize resonant quasiparticle states, so that Eq. (fl~2j) becomes: 



[u(k n ,r) 2 + v(k n ,r) 2 ] dr 



C+H+ n {k u (R- 



')) 







r + oc 




" dx + 








Jo 



CtHt (k v (R + xe i6 «)) 



dx = 1, 



(13) 



where 8 U > and 8 V > are chosen such that improper integrals converge. Hence, as in the single-particle case, 
normalization of quasiparticle S-matrix poles presents no other difficulty. As in Ref . [8[ , complex-scaled integrals will 



be denoted Reg 







f(r) dr , i.e. the regularized value of the diverging integral. 
Scattering quasiparticle states must be orthonormalized with the Dirac delta distribution: 



[u(k a , r)u(kb,r) dr + v(k a ,r)v(k b , r)] dr = 8(k a — kb), (14) 
for those with momenta k a and k\>. From Eqs. (fTC))) and (TlTj) , assuming that Eq. ^ is obtained for r > R, Eq. (I] 
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becomes: 



[u(k a , r)u(k b ,r) + v(k a ,r)v(k b ,r)] dr 



C+C+Reg 



+ C+C+Reg 



+ C Ua C^ b 



8(k a - k b ). 



C u C Uh Reg 



H tru, ( k v« r ) H tn„. (k Ub r) dr 



J R KJ k ^ r ) H LS k ^ r ) dr 

X) 

H+Jk Va r)H+ Vb (k Vb r) dr 

/+oo 
H LS^ a r)H^ Ub (k Ub r) dr 



(15) 



The divergence of the Dirac delta function at k a = k b occurs by way of the two last integrals of Eq. (fT5")) . as no 
complex scaling can make them converge if k a = k b |8| • The Dirac delta normalization of the Coulomb wave functions 
implies, as in the single-particle case: 



/-t-oo r+oo 
H ivu a ( k uar)H+ Ub (k Ub r) dr + C+ a C~ J (k Ua r)H^ (k Ub r) dr 

2irCi C u 8(k Ua — k Ub ) + f(k Ua , k Ub ), 



(16) 



where f(k Ua , k Ub ) is finite for all (k Ua , k Ub ). The relation between 5(k a — k b ) and 5(k Ua — k Ub ) is easily obtained from 
Eq. Q: 



3{k Ua k u b ) — 



dk Ua 
dk a 



(ka) 



k 

S(k a - k b ) = -j^-S(k a - k b ). 

fan. 



(17) 



This a direct application of the standard Dirac delta distribution property stating that S(f(k)) = f'(ko) 1 5(k — ko) 
for a given function f(k) bearing a unique simple zero at k = ko [291 ]. Note that k b is fixed while k a is varied to obtain 
Eq. ([T7D. Inserting Eqs. QU) and (TXTJ) into Eq. JTSJ), one obtains: 



2irk Ua 

k a 



[u(k a ,r)u(kb, r) dr + v(k a ,r)v(k bl r)] dr = 5(k a - k b ) 
Ci a C~ a 8(k a - kb) = 5(k a - kb) + g(k a , h), 



(18) 



where g{k a , k b ) bears the same properties as f(k Ua , k Ub ). As quasiparticle scattering states are orthogonal for k a ^ kb, 
g(k a , kb) = therein, so that S(k a — k b ) + g(k a , k b ) = S(k a — k b ) in all cases. 

Dirac delta distribution normalization for scattering states \k) and \k') immediately follows: 



(k\k') = S(k-k')^C+C- 



2nk„ 



(19) 



Hence, besides the additional factor k/k u , the normalization condition for quasiparticle scattering states is the same 
as that for single-particle scattering states Q. 



C. Completeness of quasiparticle states of real and complex energy 

The HFB supermatrix defined in Eq.([T]) is hermitian, so that it possesses a spectral decomposition [30l ] : 

}^ [u(k n , r)u(k n , r') + v(k n , r)v(k n ,r')] + / [u(k, r)u(k, r') + v(k, r)v(k, r')] dk = S(r — r'), (20) 

neb Jk * 

where k n — ^/2mE n /h for a bound quasisparticle state with energy E n , k is a linear momentum for a continuum 
quasiparticle state, u[k, t),v(k, r) (k — k n or k) are respectively the upper and lower components of a quasiparticle 
wave function with quantum numbers I and j (here implicit), and k\ = y— 2mX/h. All quasiparticle states must 
be normalized to one (bound) or to a Dirac delta (scattering) (see Sec. (|III Bjl ). Eq. (f20| can also be demonstrated 
extending the method of Ref. [3l| to quasi-particle states. 
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In order to obtain Berggren completeness of quasiparticle states, one can proceed as in Ref. (3^|, deforming the real 
energy contour in the complex plane. Resonant quasiparticle states appear therein, due to the Cauchy theorem, as 
S-matrix poles [32 ]. Hence, Eq. (|20|) becomes after contour deformation: 

V [u(k n ,r)u{k n ,r') + v(k„,r)v(k„,r')} + / [u(k, r)u(k, r') + v(k, r)v(k, r')] dk = 6{r - r'), (21) 

where k n refers now to a bound (b) or resonant id) (decaying) quasiparticle state and k is complex as it follows the 
deformed contour in the complex plane, denoted as L + . Resonant quasiparticle states are normalized with complex 
scaling (see Sec. piIBp ). 



IV. NUMERICAL DETERMINATION OF QUASIPARTICLE ENERGIES AND WAVE FUNCTIONS 

WITH DIRECT INTEGRATION 



A. Quasiparticle Jost functions 

In Eqs. ((9]), (|10p and (jTTJ) , constants and momenta of S-matrix poles are determined by the requirement of continuity 
of both the u(k, r) and v(k, r) functions and associated derivatives. These conditions can be expressed in a form of 
quasiparticle Jost functions, defined as a generalization of the Jost function for single-particle problems, whose zeros 
correspond to S-matrix poles (33j . They read: 



C^C+J u{k,R+) u(k,Ro) ' 
/ Cl C+\ = v'(k,R+) _ v'(k,Ry) 
V \ 'CZ'C+J v(k,R+) «(fc,i? -)' 

^V'Crci) - u(k,R^) v( kl R Q y Kll) 

where Rq is a radius typically chosen around the nuclear surface and one can demand arbitrarily that C® = = 1 
in Eqs. (HJ) and (jTUJ) . The functions, u(A;,i?o"), v(k,R^) and their derivatives, are obtained by forward integration of 
Eq. © using Eq. © as initial conditions, while .Rjj"), v(k,RQ ) and their derivatives are calculated by backward 
integration of Eq. ([2|) from the initial conditions provided by Eqs. (fTQ| and (jTTJ) . In Eq. (|22|) . one can clearly see that 
u{k,r) and v(k,r) will have continuous logarithmic derivatives if J u = and J v — respectively. However, these two 
equalities are not sufficient to uniquely determine the quasiparticle state. Indeed, they imply that one can choose 
a set of constants so that either u(k, r), u'(k, r), or v(k,r),v'(k,r) are continuous, but not necessarily both of them. 
The condition J m = is thus enforced in Eq. (|22p . The set of three equations, J u = 0, J v — and J m = 0, uniquely 
determine quasiparticle S-matrix poles. 

For quasiparticle scattering states, the linear momentum k is fixed, but constants have to be calculated with a 
matching procedure. One starts with imposing the condition = 1, as for S-matrix poles. As the u(k, r) component 
is of scattering type, the condition J u = can always be fulfilled with appropriately chosen C+ and C~ constants. 
Thus, it is sufficient to deal only with J v and J„ 



1 m • 



, ,C° V r+ \ _ v'(k,R+) v'(k,Ro) 



C°' "J „(4,<) „(*,RJ)' 
j ,% C A = M^-HMi> (23, 



,C°' v J u(k,R ) v(k,R 

the difference with Eq. (122p being that J v and J m now depend on two parameters instead of three. As in the S-matrix 
pole equations, u(k, Rq), v(k, Rq) and their derivatives are generated by forward integration of Eq. (|2|). Concerning 
the implementation of u(k,RQ), v(k,R^) and their derivatives, however, one first continues integrating forward in 
order to obtain u(k, R),u'(k, R), R being in the asymptotic region. At this point R, u(k, R),u' (k, R) provide an 
initial condition for backward integration, while Eq. (jlip is used to initialize v(k, R),v'(k, R). In this way, we obtain 
M(fc,i?^"), w(fc, i?^~) and their derivatives. Thus, the equations J v = and J m — provide the matching constants 
rendering v(k, r), t/(fc, r) continuous. 

The conditions, J u — (for S-matrix poles), J v = and J m = 0, form a system of non-linear equations of two 
or three dimensions. Consequently, it has to be solved with multi-dimensional Newton method. The only problem 
therein is to find a good starting point from where one can attain fast convergence to the exact solution in a stable 
way. 
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B. Determination of quasiparticle energy and integration constants 

Following Ref . [2(| , it is convenient to introduce linearly independent solutions of Eq. ([2]) in order to determine the 
constants denned in Eqs. ©, (TOJ and (fTTj) : 

U ^ - C° u ( fu ° ) + C° ( fv ° ) , (24) 



' M ^ C« (l U+ )+Cu ( l u ~)+Ct( { v+ ), (25) 



v J \ 9u+ J \9u- J \ 9v+ 

where the introduced basis functions verify: 

fu (r) ~ r i+1 , g V0 (r) ~ r e+1 , f Vo (r) ~ D r e + 3 , g Uo {r) ~ -lV +3 , r - 0, 

/„±(r) ~ Hf Vu {k u r) , g v +(r) ~ H+^(k v r) , f v +{r) -> , g u ± (r) -> , r -> +oo, (26) 

with 

_ m*(0)W(0) 

^° " (2£ + 3)fi2 ' (27) 

Eq. ([2T)l is obtained inserting u(r) = r^ +1 and u(r) = — D r l+3 in the second equality of Eq. @ and solving the 
equation keeping only dominant terms. 

As the basis functions of Eqs. (|24|) and (j25|) depend only on k of the quasiparticle state, they can be calculated 
with direct integration, in a forward direction for Eq. (|2"4"|) and in a backward direction for Eq. (|23|) . Used methods 
to determine quasiparticle wave function differ according to their characters; S-matrix poles or scattering states, as 
discussed below. 



C. Bound and resonant quasiparticle states 

To find S-matrix poles, it is first necessary to start with a good approximation of k, denoted k app . For that, a 
no-pairing approximation is firstly performed. Neglecting h in Eq. |T]), the Gamow-HFB equations reduce to the GHF 
equations: 

hlfa) = ei\<t>i), (28) 

where e, are complex (real) for resonant (bound) states. Eq. (|28[) provides bound and narrow resonant single-particle 
states of interest, which will be in finite number. As pairing potential h is weak compared to h, there will always 
be unique correspondence between the GHF single-particle S-matrix poles and the HFB quasiparticle S-matrix poles. 
Unless the quasiparticle S-matrix poles lie close to the Fermi energy, their lower (upper) components will be very 
close to 4>i(r) if \4>i) are (un)occupied at the HF level, so that the auxiliary energies e,, defined in Eq. 0, will be very 
close to the real parts of ej. Secondly, the HFB matrix in Eq. ([T]) is diagonalized. It has been found that the use of a 
Poschl-Teller-Ginocchio (PTG) basis provides sufficiently precise results [19| . Therefore, for Ei in Eq. ([7]) we use the 
quasiparticle energies obtained by diagonalizing the HFB matrix in the PTG basis. For a given GHF state of energy 
the starting quasiparticle energy E app (from which k app is immediately deduced), is then the BCS quasiparticle 
energy whose e~j is closest to the real part of e; . If the HFB quasiparticle S-matrix pole is far from the Fermi energy, 
E a pp is very close to the exact energy. Otherwise, it will still provide a good starting point, as, in practice, one can 
have only one quasiparticle state close to the Fermi energy for a given (£, j)-partial wave. 

Furthermore, one demands C~ = 0, which translates into a linear eigen-value problem of dimension equal to four, 
deduced from Eqs. (j2"4")l and ([23]) . which one matches at r — Rq: 

- 0, (29) 

where all matrix functions have been evaluated at r = Rq by way of backward or forward integration. As the 
integration constants are not simultaneously equal to zero, they have to form an eigenvector of the matching matrix 
of Eq. (|2"9"|) , which we denote M hereafter, of eigenvalue equal to zero. However, the determinant of the 4x4 matrix 
M is zero uniquely for the exact value of k. Thus, the set of approximate constants to use as a starting point for 
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Newton method is chosen as the eigenvector of l MM whose associated eigenvalue is the smallest in modulus ( l MM 
is used instead of M because it is symmetric). The constant ratios C°/C° and C+/C+ used in Eq. (|2"2")l follow, as 
they are independent of the norm of the considered eigenvector. Exact determination of k, C°/C° and C+/C+ can 
then be worked out via three-dimensional Newton method. 



D. Scattering quasiparticle state 

If one considers a scattering state, it is convenient to define a + , a" , b + , b~ so that — a ± C^ l + b C°. Moreover, 
as all constants are calculated up to a normalization factor, one can impose C° = 1. Upper components of Eqs. (|24|) 
and (|25[) matched at r = R and Eq. (j2l)|) provide linear equations for a ± and 6 ± : 

a+f u+ (R)+a-f u -(R) = f ua (R) , b+ f u+ (R) + b~ f a - (R) = f V0 (R), 

a+f' u+ (R) + a-f' u - (R) = f Uo (R) , b+ f' u+ (R) + b~ f' u - (R) = f' Vo (R). (30) 

From the knowledge of a* and b ± , matching lower components in Eqs. (|24|) and |25|) at r = Ro determines C° and 
C+ via linear equations as well: 

C v[dv ( R o) - b + g u +(R ) - b~g u -(Ro)] - C+g v +{R ) = a+ g u +(R Q ) + a~ g u -(R ) - g Uo (Ro), 

C° v K ( R o) - b + g' u+ (Ro) - b-g' u -(Ro)} - C+g' v+ (Ro) = a+g' u+ (Ro) + a -g' u -(R ) - g' Uo (R ). (31) 

As — a ± + 6 ± C°, all constants are determined with simple two-dimensional linear systems. Newton method 
applied to Eq. ([2"3"]l converges very quickly using the obtained set of constants as a starting point. Note that the use 
of Hf v {k u r) functions in Eqs. (|2"5|) and (f2T))) can be sometimes unstable, especially for the proton case, where, for 
low scattering energies, they can be very large and cancel almost exactly in Eq. (flUl) . In this case, it is better to use 
regular and irregular Coulomb wave functions, Fi Vu (k u r) and Gi riu (k u r), as basis functions. 



V. NORMAL AND PAIRING DENSITIES 



As quasiparticle states of complex energy form a complete set (see Eq. (|2~T]) ), one can directly express densities with 
upper and lower components of quasiparticle states: 

n£(M) t-3 

Mr) = - £ - /*r W r) * , «,) = ^fcM. W 

ne(6,d) i+ tj 

where Pij(r) and f>ij(r) are respectively partial normal and pairing densities related to a given partial wave with 
quantum numbers i and j, and p(r), p(r) are respectively the normal and pairing densities of the HFB ground state. 
However, due to the zero-range character of Skyrme forces, it is necessary to introduce an energy cut in contour 
integrals, so that L + contour has to stop at finite energy E cut (see Fig. [T]). Note that, due to this requirement, 
it is necessary for L + complex contours to come back to the real axis. Even though quasiparticle wave functions 
are complex in Eq. (|32p , pij(r) and pijir) are real because one is considering a HFB bound ground state, so that, 
due to Cauchy theorem, complex integration in Eq. (|32[) is equivalent to real integration in the standard case. As a 
consequence, the DFT can be applied also to the Gamow HFB formalism, i.e. potentials V(r) and W(r) in Eq. are 
evaluated using the standard formulas of Ref. [ijj]. As shown in Fig. [I] the bound HF single-particle states can become 
resonant states when pairing correlations are added [37| . Thus, physical interpretation of a resonant quasiparticle is 
somewhat different from that of single-particle resonances: widths of the quasiparticle states associated with the HF 
bound single-particle states originates from pairing-induced couplings between the bound and scattering states [3~jj . 

In the same way as in the Gamow Shell Model [JH, 0] , the scattering L + contours in Eq. (|3^|) have to be discretized, 
providing a finite set of linear momenta and weights (ki,Wi). In practice, the Gauss-Legendre quadrature has been 
found to be most efficient. Scattering quasiparticle states are also renormalized, multiplying them by y/wi [32| . so 
that the discretized expressions of Eq. (j32|) are formally identical to the discrete case: 

n£(b,d) i 

pij{r) ~ - 2_j u(k ni r)v(k n ,r) ~^2u Wi {ki,r)v Wi (ki,r), (33) 

ne(b,d) i 
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where u Wi (k i: r) = ^Jwl it(fe;,r) and v Wi (k i: r) = ^/wl v(ki,r). 



VI. ANOTHER METHOD: EXPANSION OF QUASIPARTICLE STATES WITH THE GHF BASIS 



Another possibility to solve the HFB equations in complex energy plane is to use the Gamow single-particle states 
as a basis. The optimal Berggren basis to expand the HFB quasiparticle states is obviously the GHF basis generated 
by the potential V(r) and the effective mass m*{r) of Eq. @. Note that it is not equivalent to the GHF basis issued 
from the pure HF variational principle in that pairing correlations always give extra contributions to the particle-hole 
part of the HFB Hamiltonian. Indeed, we noticed in our numerical calculation that other Berggren bases make the 
HFB self-consistent procedure unstable due to the appearance of very large matrix elements in the HFB Hamiltonian 
matrix. The use of the optimized Berggren basis mentioned above removes this problem. This approach may be 
regarded as a generalization of the two-basis method [HI, H3, [Hj] . 

The GHF basis states <p(r) are defined by the following equation: 



dr 2m* (r) dr 



n 2 e(e + i) 

2m*(r)r 2 



V(r) 



')> 



(34) 



issued directly from Eqs. and ([28]). where e is the complex energy of the GHF state. The HFB Hamiltonian matrix 
represented with this basis becomes: 



h- X h 
h X-h 



ex - X 


\ 




h 


ejy — A 






A-ei 


h 






A - ejv / 



(35) 



V 

where the continuous Berggren basis is discretized with the Gauss-Legendre quadrature (see Sec. |VJ) so that total 
number of basis states is N. Its particle-hole part is evidently diagonal and matrix elements of h read: 



r+oo 

(<t>a\h\4>b) = / Mr)W(r)Mr) dr, (36) 
Jo 

where \4>b) and \cj> a ) are the GHF basis states and W(r) is the HFB pairing potential defined in Eq. (|2|). For bound 
HFB ground states, W(r) decreases sufficiently quickly so that no complex scaling is needed to evaluate the integral of 
Eq. (|3"6"1) . Hence, after discretization of the contours representing scattering basis states, this method takes a formally 
identical form to the standard matrix diagonalization treatment of the HFB problem. 



VII. NUMERICAL APPLICATIONS 



The frameworks described above, i.e. the Gamow-HFB approach in the coordinate or the GHF configurational 
space, are applied to Nickel isotopes close to the neutron drip-line, from 84 Ni to 90 Ni, which possess spherical HFB 
ground states. In the numerical calculation, the SLy4 Skyrme force [35} is used in combination with the surface-type 
contact pairing interaction [26j] whose pairing strength is fitted to reproduce the pairing gap of 120 Sn. Using the 
standard notation [2||, the pairing interaction parameters read t = —519.9 MeV fm 3 for the density-independent 
part and t 3 — — 37.5i MeV fm 6 for the density-dependent part. The maximal angular momentum used is l max = 10 
and a sharp cut-off at E cut = 60 MeV is adopted. Scattering contours of quasiparticle states are discretized with 60, 
100 or 300 Gaussian points. Several hundred points are indeed necessary when resonant states lie relatively close to 
E cut (see Fig. [T]), as is the case for the HFB quasiparticle resonance associated with the deeply bound neutron 0si/ 2 
HF state for example (see Table H| . Scattering contours of single-particle states in the GHF basis are discretized up 
to k max = 4 fm -1 with 100 points, which in this case assures convergence of numerical calculation. This concerns 
only for the neutron channel, as the pairing gap vanishes in the proton channel. 

The result of calculation for normal and pairing densities are presented in Figs. [Hand [3] It is interesting to compare 
the densities obtained by solving the Gamow-HFB equations in the coordinate or the GHF configurational space to 
those calculated by the standard coordinate space framework where the continuum is discretized with box boundary 
conditions. They are denoted GHFB/Coord., GHFB/Config. and HFB/Box., respectively. All results coincide in 
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both normal and logarithmic scales for r < 30 fm. It was also checked that spurious imaginary parts of densities, 
caused by the discretization of the continuum of complex energy, were negligible, of the order of 10 -6 [fm -3 ] for 
GHFB/Coord. and 10 -12 [fm -3 ] for GHFB/Config., as the largest error values. In Tabled the bound and resonant 
single-particle states obtained by the GHF calculation are compared with the corresponding quasiparticle states 
calculated by the GHFB/Coord. method. It is obvious that bound HF states can give rise to unbound quasiparticle 
states carrying a sizable width when pairing correlations are switched on. 

Physical observables associated with the HFB ground states are provided in Tables [Til an d IIIII On the one hand, 
differences occur for neutron pairing energies, which are most sensitive to continuum effects [36j . While those of 
GHFB/Coord. compared to HFB/Box remain of the order of 500 keV, the difference between GHFB/Config. and 
HFB/Box pairing energies can be ~1.5 MeV. On the other hand, the r.m.s. radii and total energies are basically 
the same, with a discrepancy of at most ~300 keV for the latter. These results indicate that the GHFB/Coord., 
GHFB/Config. and HFB/Box treatments are all reliable methods to solve the HFB equations taking the continuum 
effects into account. As resonant states are explicitly treated in the Gamow HFB approach, this implies that the 
resonant effects can be well accounted for also by means of the HFB/Box method. This point is not necessarily 
widely accepted [13] ■ Even though the good agreements among the results of the GHFB/Coord., GHFB/Config. and 
HFB/Box calculations might be surprising, we see no reason to suspect that this is an exceptional case valid only for 
the Ni isotopes considered here. It will be interesting to examine this point further. 

VIII. PERSPECTIVES FOR DESCRIBING DECAYING NUCLEI AND BEYOND-MEAN FIELD 

APPROACHES 

The GHFB/Coord. method directly provides quasiparticle wave functions without using any intermediate basis 
states. Hence, it may be used also to describe decaying nuclear ground states in the HFB approximation. In fact, 
no HFB theory capable of describing decaying HFB ground states exists, even though an approximate scheme was 
proposed in Ref. [38|]. The main difficulty is that it is not possible to construct the HFB ground state obeying the 
outgoing wave condition if one includes the full set of quasiparticle states of positive energy 38]. This arises from 
the fact that quasiparticles form a degenerate continuum of scattering states for E < |A| if A > 0, whereas they can 
only generate a discrete set of bound states in this region if A < 0. It is impossible to remove quasiparticle states 
with E < |A| with the use of the GHFB/Config. method, because quasiparticle eigen-energies of the HFB matrix are 
complex. In contrast, the direct integration method (GHFB/Coord.) allows us to select which quasiparticle states 
are occupied in the HFB ground state. Hence, it may be possible to carefully study properties of decaying HFB states 
at least for the spherical case. 

The GHF configurational approach (GHFB/Config.) may be more appropriate to study excited states in deformed 
nuclei by means of the QRPA. For deformed nuclei, basis expansion approaches may be easier compared to the 
calculation of deformed HFB ground states in coordinate space [l(|. For calculating bound HFB ground states, we 
can use the PTG basis, which is more efficient than the GHF basis, considering the numerical cost of recalculating the 
GHF basis states inherent to the two-basis method (see Sec. IVI|) . Once a HFB ground state is obtained in this way, 
one can readily calculate the GHF basis wave functions. The QRPA matrix would then be represented afterward with 
respect to the quasiparticles wave functions expanded in the GHF basis, thus allowing the description of unbound 
QRPA excited states. 

IX. CONCLUSION 

The Berggren completeness relation, originally developed in the context of standard Schrodinger equation, has been 
extended to quasiparticles in the HFB formalism. It was shown that, as in the standard single-particle potential 
problem, bound, resonant and scattering quasiparticles are well defined and form a complete set, by which bound 
HFB ground states can be constructed. Both situations are very similar and can be treated by contour deformation of 
continuous real sets of states, even though physical interpretation of resonant quasiparticles is different from that of 
resonant single-particles. Numerical applications have been effected with neutron- rich Nickel isotopes close to the drip 
line, for which continuum coupling is important. It was shown that the Gamow-HFB approach, both in coordinate 
and configurational space representations, properly describe densities and physical observables. Thus, it provides us 
with an efficient tool to study ground states of medium and heavy nuclei close to the drip line. With these approaches, 
QRPA calculation fully taking into account continuum coupling may be efficiently carried out. 
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TABLE I: Bound and resonant neutron energies and widths for Ni, calculated in the GHF approximation and in the 
GHFB/Coord. formalism. Single-particle energies (ej) and quasiparticle energies (Ei) are given in MeV and widths (T) in 
keV. Note that the GHF 2si/ 2 state dissolves into continuum quasiparticle states in the Gamow-HFB description. 



GHF GHFB/Coord. 



states e T E V 



OSl/2 




-52.618 





51.573 




1.099 10 


lSl/2 




-24.630 





24.348 




46.006 


2si/ 2 




-1.196 











0P3/2 




-41.655 





40.796 




27.282 


1P3/2 




-12.986 





12.658 




490.565 


Opi/2 




-42.881 





38.870 




27.138 


lPl/2 




-11.189 





10.816 




404.299 


0d 5/2 




-29.921 





29.141 




0.780 


la 5/2 









3.181 




194.181 


U"3/2 




-30.657 





25.095 




22.567 


1^3/2 




-0.349 





2.173 




560.608 


0/7/2 




-18.177 





17.654 




397.374 


0/5/2 




-11.331 





11.065 




645.638 


039/2 




-6.770 





6.570 




0.807 


037/2 




1.350 


6.410 


3.120 




63.6131 


0/lll/2 




3.852 


52.851 


5.269 




131.776 


TABLE 


II: Gamow-HFB 


■ observables for 


84 Ni and 86 Ni calculated with the GHFB/Coord., GHFB/Config. and HFB/Box 


methods. 


The r.m.s. radii are given in fm and other quantities 


in MeV. The proton 


chemical potential X p is not presented as 


there is no proton pairing gap. 














84 Ni 






86 Ni 






HFB/Box 


GHFB/Coord. 


GHFB/Config. 


HFB/Box GHFB/Coord. 


GHFB/Config. 


An 


-1.453 


-1.430 


-1.440 


-1.037 


-1.027 


-1.029 


r n 


4.451 


4.450 


4.450 


4.528 


4.526 


4.526 


r p 


3.980 


3.982 


3.982 


4.001 


4.001 


4.001 


A n 


1.481 


1.535 


1.564 


1.667 


1.658 


1.669 


E pair 


-30.70 


-30.72 


-31.85 


-36.52 


-35.85 


-36.39 


T n 


1084.53 


1086.05 


1086.46 


1118.65 


1118.68 


1118.78 




430.47 


430.23 


430.17 


425.99 


426.01 


426.00 


TpSO 


-63.379 


-63.164 


-63.01 


-61.679 


-61.712 


-61.631 


TpCoul 


132.94 


132.89 


132.88 


132.26 


132.25 


132.25 


^exc 


-10.138 


-10.135 


-10.135 


-10.084 


-10.085 


-10.085 


Etot 


-654.89 


-654.89 


-655.05 


-656.933 


-656.836 


-656.971 
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TABLE III: Same as in Table El but for S8 Ni and 90 Ni. 

88 Ni 90 Ni 
HFB/Box GHFB/Coord. GHFB/Config. HFB/Box GHFB/Coord. GHFB/Config. 



An 


-0.671 


-0.661 


-0.665 


-0.342 


-0.330 


-0.342 


r„ 


4.603 


4.602 


4.601 


4.677 


4.674 


4.675 


r p 


4.021 


4.022 


4.022 


4.043 


4.043 


4.043 


A n 


1.790 


1.782 


1.800 


1.899 


1.899 


1.935 




-41.98 


-41.26 


-42.17 


-47.158 


-46.509 


-48.449 


T„ 


1150.71 


1150.74 


1151.02 


1182.52 


1182.91 


1183.79 




421.71 


421.71 


421.70 


417.38 


417.35 


417.31 




-59.558 


-59.559 


-59.470 


-56.898 


-56.887 


-56.822 


rpCoul 

' dir 


131.571 


131.576 


131.576 


130.947 


130.883 


130.878 


^eic 


-10.033 


-10.033 


-10.033 


-9.980 


-9.980 


-9.980 


Etot 


-658.167 


-658.082 


-658.272 


-658.665 


-658.635 


-658.936 



Im(k) 



bound quasi-particle 



a |A,| 



1/2 



resonant quasi-particle 
(HF bound or resonant) 



\ 



cut 



Re(k) 



L contour 

(scattering states) 



FIG. 1: (color online) Location of quasiparticle S-matrix poles and deformed complex contour L + of scattering quasiparticle 
states used in the Berggren completeness relation. Here, a — \/2m/h. 
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r (fm) r (fm) 



FIG. 2: (color online) Neutron densities p n both in normal (left-hand side) and logarithmic (rightdiand side) scales. Results of 
the HFB/Box, GHFB/Coord. and GHFB/Config. calculations are displayed by solid, dashed and dotted lines, respectively. 



